week 2: linear model and causal inference

geocentric models

Annoucements and such

  • We will start using brms today! Install this package now, if you haven’t already.
install.packages(c("brms","tidybayes"))

Workspace setup

library(here)
library(tidyverse)
library(cowplot)
library(brms)
library(tidybayes)
library(patchwork)

Workflow

  1. State a clear question.

  2. Sketch your causal assumptions.

  3. Use the sketch to define a generative model.

  4. Use the generative model to build an estimator.

  5. Profit.

Model “recipes”

  1. Recognize a set of variables to work with. (Data and parameters.)
  2. Define each variable either in terms of the other variables OR in terms of a probability distribution.
  3. The combination of variables and their probability distributions defines a joint generative model that can be used to simulate hypothetical observations and analyze real ones.

Here’s an example:

\[\begin{align*} y_i &\sim \text{Normal}(\mu_i,\sigma) \\ \mu_i &= \beta x_i \\ \beta &\sim \text{Normal}(0,10) \\ \sigma &\sim \text{Exponential}(1) \\ x_i &\sim \text{Normal}(0,1) \\ \end{align*}\]

Model for globe-tossing

Here’s the model for last week’s globe-tossing experiment:

\[\begin{align*} W &\sim \text{Binomial}(N,p) \\ p &\sim \text{Uniform}(0,1) \\ \end{align*}\]

  • \(W\) is the observed count of water.
  • \(N\) is the total number of tosses.
  • \(p\) is the proportion of water on the globe.

The whole model can be read as:

The count \(W\) is distributed binomially with sample size \(N\) and probability \(p\). The prior for \(p\) is assumed to be uniform between 0 and 1.

Model for globe-tossing

Here’s the model for last week’s globe-tossing experiment:

\[\begin{align*} W &\sim \text{Binomial}(N,p) \\ p &\sim \text{Uniform}(0,1) \\ \end{align*}\]

Estimating the posterior using brms

Last week, we used grid approximation to estimate the posterior distribution. In the video you watched for today, McElreath moves on to something called QUADRATIC APPROXIMATION. It’s good to understand what that’s doing, but you and I are moving right along to MCMC. We won’t go into the details of what’s happening for a few weeks, but let’s start with the code to estimate the model.

Let’s say we tossed the globe 9 times and observed 6 waters:

m1 <-
  brm(data = list(w = 6),
      family = binomial(link = "identity"),
      w | trials(9) ~ 0 + Intercept,
      prior(uniform(0, 1), class = b),
      iter = 5000, warmup = 1000, seed = 3, chains=1,
      file = here("files/models/m21.1"))
1
Data can be a data frame or a list. Just make sure the variable names match what’s in your formula.
2
How you assume your outcome variable is distributed.
3
The formula for your outcome.
4
Priors for every parameter in your model. Here, we only have one parameter, so we only need 1 prior. This happens to be a flat prior between 0 and 1.
5
Some choices about how we want our model to run. We’ll go more into this later.
6
These models can take a long time to run. You can automatically save the output to a file; when you do this, the next time you run this code, it won’t actually estimate the model, but will instead pull the output from your stated file. Be WARNED: if you change the data or the model code, it will NOT restimate your model until you delete the file.

Sampling from the posterior

Grid approximation gave us the calculated probability of each possible value of our parameter, \(p\). But our method of conducting bayes will no longer give us such a neat solution. Here’s how you get the posterior distribution for \(p\):

samples_from_post = as_draws_df(m1)
samples_from_post
# A draws_df: 4000 iterations, 1 chains, and 3 variables
   b_Intercept lprior lp__
1         0.46      0 -2.1
2         0.49      0 -1.9
3         0.56      0 -1.5
4         0.66      0 -1.3
5         0.71      0 -1.3
6         0.62      0 -1.3
7         0.68      0 -1.3
8         0.55      0 -1.5
9         0.63      0 -1.3
10        0.67      0 -1.3
# ... with 3990 more draws
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}
samples_from_post %>%  
  ggplot(aes(x=b_Intercept)) +
  geom_density(fill = "grey", color = "white") +
  labs(x="Proportion water")

Posterior predictive distribution

ppd = posterior_predict(m1)
dim(ppd)
[1] 4000    1
ppd
        [,1]
   [1,]    5
   [2,]    5
   [3,]    6
   [4,]    5
   [5,]    6
   [6,]    4
   [7,]    8
   [8,]    4
   [9,]    7
  [10,]    8
  [11,]    6
  [12,]    5
  [13,]    5
  [14,]    8
  [15,]    8
  [16,]    7
  [17,]    7
  [18,]    9
  [19,]    7
  [20,]    8
  [21,]    8
  [22,]    8
  [23,]    5
  [24,]    5
  [25,]    7
  [26,]    6
  [27,]    6
  [28,]    5
  [29,]    5
  [30,]    5
  [31,]    4
  [32,]    7
  [33,]    8
  [34,]    6
  [35,]    3
  [36,]    4
  [37,]    3
  [38,]    7
  [39,]    6
  [40,]    8
  [41,]    6
  [42,]    7
  [43,]    8
  [44,]    4
  [45,]    6
  [46,]    5
  [47,]    7
  [48,]    4
  [49,]    8
  [50,]    9
  [51,]    5
  [52,]    4
  [53,]    4
  [54,]    6
  [55,]    5
  [56,]    7
  [57,]    4
  [58,]    0
  [59,]    7
  [60,]    6
  [61,]    6
  [62,]    7
  [63,]    8
  [64,]    3
  [65,]    7
  [66,]    5
  [67,]    9
  [68,]    3
  [69,]    6
  [70,]    8
  [71,]    6
  [72,]    3
  [73,]    5
  [74,]    7
  [75,]    6
  [76,]    3
  [77,]    8
  [78,]    7
  [79,]    8
  [80,]    4
  [81,]    4
  [82,]    5
  [83,]    7
  [84,]    8
  [85,]    6
  [86,]    3
  [87,]    3
  [88,]    4
  [89,]    7
  [90,]    7
  [91,]    8
  [92,]    5
  [93,]    5
  [94,]    3
  [95,]    4
  [96,]    9
  [97,]    9
  [98,]    6
  [99,]    9
 [100,]    8
 [101,]    4
 [102,]    7
 [103,]    7
 [104,]    6
 [105,]    4
 [106,]    7
 [107,]    7
 [108,]    7
 [109,]    6
 [110,]    7
 [111,]    6
 [112,]    4
 [113,]    5
 [114,]    6
 [115,]    6
 [116,]    3
 [117,]    7
 [118,]    4
 [119,]    7
 [120,]    5
 [121,]    6
 [122,]    5
 [123,]    8
 [124,]    8
 [125,]    8
 [126,]    6
 [127,]    8
 [128,]    6
 [129,]    6
 [130,]    6
 [131,]    8
 [132,]    5
 [133,]    4
 [134,]    8
 [135,]    5
 [136,]    5
 [137,]    6
 [138,]    6
 [139,]    8
 [140,]    9
 [141,]    5
 [142,]    4
 [143,]    7
 [144,]    5
 [145,]    5
 [146,]    7
 [147,]    5
 [148,]    4
 [149,]    7
 [150,]    5
 [151,]    4
 [152,]    8
 [153,]    7
 [154,]    5
 [155,]    7
 [156,]    6
 [157,]    5
 [158,]    7
 [159,]    7
 [160,]    4
 [161,]    8
 [162,]    6
 [163,]    4
 [164,]    3
 [165,]    3
 [166,]    6
 [167,]    3
 [168,]    4
 [169,]    6
 [170,]    7
 [171,]    4
 [172,]    1
 [173,]    6
 [174,]    9
 [175,]    5
 [176,]    3
 [177,]    9
 [178,]    6
 [179,]    8
 [180,]    8
 [181,]    3
 [182,]    8
 [183,]    7
 [184,]    7
 [185,]    7
 [186,]    6
 [187,]    7
 [188,]    8
 [189,]    2
 [190,]    4
 [191,]    7
 [192,]    7
 [193,]    8
 [194,]    4
 [195,]    5
 [196,]    6
 [197,]    6
 [198,]    7
 [199,]    6
 [200,]    8
 [201,]    6
 [202,]    8
 [203,]    9
 [204,]    9
 [205,]    3
 [206,]    7
 [207,]    6
 [208,]    4
 [209,]    9
 [210,]    5
 [211,]    9
 [212,]    3
 [213,]    8
 [214,]    7
 [215,]    7
 [216,]    7
 [217,]    8
 [218,]    7
 [219,]    6
 [220,]    7
 [221,]    7
 [222,]    7
 [223,]    3
 [224,]    9
 [225,]    7
 [226,]    6
 [227,]    3
 [228,]    8
 [229,]    7
 [230,]    8
 [231,]    7
 [232,]    4
 [233,]    7
 [234,]    9
 [235,]    6
 [236,]    7
 [237,]    5
 [238,]    7
 [239,]    5
 [240,]    3
 [241,]    4
 [242,]    8
 [243,]    6
 [244,]    5
 [245,]    5
 [246,]    7
 [247,]    3
 [248,]    8
 [249,]    5
 [250,]    3
 [251,]    9
 [252,]    7
 [253,]    7
 [254,]    2
 [255,]    5
 [256,]    8
 [257,]    5
 [258,]    5
 [259,]    6
 [260,]    4
 [261,]    7
 [262,]    5
 [263,]    4
 [264,]    5
 [265,]    4
 [266,]    5
 [267,]    5
 [268,]    5
 [269,]    6
 [270,]    6
 [271,]    5
 [272,]    7
 [273,]    5
 [274,]    5
 [275,]    8
 [276,]    4
 [277,]    5
 [278,]    8
 [279,]    8
 [280,]    7
 [281,]    5
 [282,]    6
 [283,]    3
 [284,]    5
 [285,]    3
 [286,]    4
 [287,]    5
 [288,]    7
 [289,]    8
 [290,]    9
 [291,]    6
 [292,]    6
 [293,]    8
 [294,]    6
 [295,]    6
 [296,]    5
 [297,]    7
 [298,]    6
 [299,]    5
 [300,]    4
 [301,]    5
 [302,]    6
 [303,]    6
 [304,]    5
 [305,]    7
 [306,]    9
 [307,]    8
 [308,]    7
 [309,]    7
 [310,]    7
 [311,]    6
 [312,]    9
 [313,]    6
 [314,]    9
 [315,]    5
 [316,]    8
 [317,]    8
 [318,]    6
 [319,]    6
 [320,]    6
 [321,]    7
 [322,]    4
 [323,]    4
 [324,]    6
 [325,]    6
 [326,]    8
 [327,]    6
 [328,]    8
 [329,]    5
 [330,]    5
 [331,]    5
 [332,]    9
 [333,]    6
 [334,]    7
 [335,]    7
 [336,]    7
 [337,]    9
 [338,]    2
 [339,]    3
 [340,]    1
 [341,]    2
 [342,]    5
 [343,]    6
 [344,]    7
 [345,]    4
 [346,]    6
 [347,]    5
 [348,]    8
 [349,]    7
 [350,]    6
 [351,]    9
 [352,]    7
 [353,]    5
 [354,]    7
 [355,]    6
 [356,]    6
 [357,]    5
 [358,]    5
 [359,]    2
 [360,]    2
 [361,]    6
 [362,]    5
 [363,]    4
 [364,]    7
 [365,]    7
 [366,]    7
 [367,]    4
 [368,]    2
 [369,]    8
 [370,]    9
 [371,]    5
 [372,]    6
 [373,]    5
 [374,]    7
 [375,]    4
 [376,]    2
 [377,]    3
 [378,]    6
 [379,]    4
 [380,]    7
 [381,]    7
 [382,]    8
 [383,]    8
 [384,]    7
 [385,]    8
 [386,]    8
 [387,]    6
 [388,]    8
 [389,]    4
 [390,]    6
 [391,]    8
 [392,]    7
 [393,]    5
 [394,]    6
 [395,]    5
 [396,]    6
 [397,]    8
 [398,]    6
 [399,]    9
 [400,]    8
 [401,]    7
 [402,]    5
 [403,]    8
 [404,]    8
 [405,]    5
 [406,]    7
 [407,]    7
 [408,]    9
 [409,]    6
 [410,]    9
 [411,]    7
 [412,]    9
 [413,]    6
 [414,]    4
 [415,]    9
 [416,]    7
 [417,]    6
 [418,]    4
 [419,]    8
 [420,]    6
 [421,]    7
 [422,]    7
 [423,]    7
 [424,]    7
 [425,]    9
 [426,]    2
 [427,]    3
 [428,]    8
 [429,]    7
 [430,]    7
 [431,]    4
 [432,]    5
 [433,]    3
 [434,]    6
 [435,]    8
 [436,]    7
 [437,]    6
 [438,]    4
 [439,]    6
 [440,]    8
 [441,]    7
 [442,]    6
 [443,]    8
 [444,]    8
 [445,]    7
 [446,]    4
 [447,]    5
 [448,]    3
 [449,]    2
 [450,]    4
 [451,]    6
 [452,]    9
 [453,]    6
 [454,]    6
 [455,]    5
 [456,]    4
 [457,]    7
 [458,]    7
 [459,]    8
 [460,]    5
 [461,]    7
 [462,]    8
 [463,]    7
 [464,]    7
 [465,]    7
 [466,]    1
 [467,]    6
 [468,]    4
 [469,]    7
 [470,]    6
 [471,]    4
 [472,]    8
 [473,]    5
 [474,]    4
 [475,]    4
 [476,]    7
 [477,]    5
 [478,]    7
 [479,]    7
 [480,]    6
 [481,]    5
 [482,]    8
 [483,]    6
 [484,]    6
 [485,]    3
 [486,]    7
 [487,]    4
 [488,]    3
 [489,]    5
 [490,]    5
 [491,]    7
 [492,]    6
 [493,]    6
 [494,]    6
 [495,]    5
 [496,]    7
 [497,]    5
 [498,]    5
 [499,]    4
 [500,]    5
 [501,]    4
 [502,]    6
 [503,]    5
 [504,]    6
 [505,]    6
 [506,]    8
 [507,]    4
 [508,]    6
 [509,]    5
 [510,]    5
 [511,]    7
 [512,]    9
 [513,]    4
 [514,]    8
 [515,]    8
 [516,]    5
 [517,]    6
 [518,]    7
 [519,]    9
 [520,]    3
 [521,]    5
 [522,]    1
 [523,]    2
 [524,]    3
 [525,]    6
 [526,]    9
 [527,]    9
 [528,]    5
 [529,]    5
 [530,]    4
 [531,]    9
 [532,]    9
 [533,]    7
 [534,]    8
 [535,]    5
 [536,]    6
 [537,]    7
 [538,]    5
 [539,]    6
 [540,]    6
 [541,]    4
 [542,]    5
 [543,]    6
 [544,]    7
 [545,]    4
 [546,]    6
 [547,]    7
 [548,]    5
 [549,]    6
 [550,]    4
 [551,]    8
 [552,]    8
 [553,]    7
 [554,]    7
 [555,]    6
 [556,]    5
 [557,]    3
 [558,]    4
 [559,]    6
 [560,]    6
 [561,]    6
 [562,]    6
 [563,]    6
 [564,]    5
 [565,]    3
 [566,]    8
 [567,]    6
 [568,]    5
 [569,]    6
 [570,]    4
 [571,]    6
 [572,]    7
 [573,]    3
 [574,]    3
 [575,]    3
 [576,]    7
 [577,]    8
 [578,]    8
 [579,]    5
 [580,]    5
 [581,]    7
 [582,]    6
 [583,]    6
 [584,]    4
 [585,]    4
 [586,]    3
 [587,]    2
 [588,]    4
 [589,]    2
 [590,]    1
 [591,]    8
 [592,]    1
 [593,]    4
 [594,]    8
 [595,]    4
 [596,]    8
 [597,]    7
 [598,]    4
 [599,]    6
 [600,]    7
 [601,]    4
 [602,]    5
 [603,]    5
 [604,]    6
 [605,]    8
 [606,]    6
 [607,]    7
 [608,]    6
 [609,]    3
 [610,]    8
 [611,]    3
 [612,]    7
 [613,]    6
 [614,]    5
 [615,]    8
 [616,]    5
 [617,]    6
 [618,]    7
 [619,]    6
 [620,]    6
 [621,]    2
 [622,]    6
 [623,]    7
 [624,]    7
 [625,]    3
 [626,]    6
 [627,]    7
 [628,]    5
 [629,]    8
 [630,]    3
 [631,]    4
 [632,]    5
 [633,]    4
 [634,]    7
 [635,]    7
 [636,]    6
 [637,]    7
 [638,]    4
 [639,]    3
 [640,]    6
 [641,]    6
 [642,]    7
 [643,]    7
 [644,]    6
 [645,]    5
 [646,]    6
 [647,]    1
 [648,]    3
 [649,]    3
 [650,]    7
 [651,]    6
 [652,]    7
 [653,]    7
 [654,]    8
 [655,]    6
 [656,]    8
 [657,]    7
 [658,]    8
 [659,]    6
 [660,]    6
 [661,]    5
 [662,]    3
 [663,]    3
 [664,]    4
 [665,]    8
 [666,]    7
 [667,]    8
 [668,]    5
 [669,]    7
 [670,]    3
 [671,]    7
 [672,]    5
 [673,]    7
 [674,]    5
 [675,]    6
 [676,]    4
 [677,]    7
 [678,]    6
 [679,]    8
 [680,]    4
 [681,]    5
 [682,]    6
 [683,]    6
 [684,]    6
 [685,]    6
 [686,]    7
 [687,]    6
 [688,]    6
 [689,]    7
 [690,]    5
 [691,]    9
 [692,]    4
 [693,]    5
 [694,]    5
 [695,]    5
 [696,]    7
 [697,]    5
 [698,]    3
 [699,]    2
 [700,]    2
 [701,]    5
 [702,]    4
 [703,]    6
 [704,]    8
 [705,]    5
 [706,]    6
 [707,]    4
 [708,]    4
 [709,]    7
 [710,]    6
 [711,]    6
 [712,]    8
 [713,]    7
 [714,]    5
 [715,]    6
 [716,]    7
 [717,]    8
 [718,]    3
 [719,]    6
 [720,]    5
 [721,]    3
 [722,]    4
 [723,]    4
 [724,]    6
 [725,]    8
 [726,]    4
 [727,]    4
 [728,]    6
 [729,]    7
 [730,]    2
 [731,]    4
 [732,]    8
 [733,]    5
 [734,]    5
 [735,]    7
 [736,]    6
 [737,]    6
 [738,]    7
 [739,]    6
 [740,]    6
 [741,]    6
 [742,]    9
 [743,]    7
 [744,]    7
 [745,]    8
 [746,]    6
 [747,]    6
 [748,]    3
 [749,]    7
 [750,]    2
 [751,]    7
 [752,]    6
 [753,]    4
 [754,]    5
 [755,]    6
 [756,]    4
 [757,]    3
 [758,]    6
 [759,]    5
 [760,]    6
 [761,]    6
 [762,]    4
 [763,]    6
 [764,]    8
 [765,]    5
 [766,]    7
 [767,]    6
 [768,]    7
 [769,]    7
 [770,]    5
 [771,]    6
 [772,]    6
 [773,]    3
 [774,]    8
 [775,]    7
 [776,]    8
 [777,]    5
 [778,]    8
 [779,]    4
 [780,]    5
 [781,]    4
 [782,]    4
 [783,]    4
 [784,]    5
 [785,]    3
 [786,]    3
 [787,]    3
 [788,]    4
 [789,]    5
 [790,]    8
 [791,]    9
 [792,]    8
 [793,]    8
 [794,]    4
 [795,]    5
 [796,]    3
 [797,]    5
 [798,]    6
 [799,]    7
 [800,]    6
 [801,]    8
 [802,]    9
 [803,]    7
 [804,]    5
 [805,]    5
 [806,]    6
 [807,]    6
 [808,]    8
 [809,]    6
 [810,]    5
 [811,]    5
 [812,]    5
 [813,]    6
 [814,]    9
 [815,]    9
 [816,]    6
 [817,]    4
 [818,]    6
 [819,]    5
 [820,]    8
 [821,]    4
 [822,]    4
 [823,]    8
 [824,]    7
 [825,]    4
 [826,]    4
 [827,]    5
 [828,]    5
 [829,]    2
 [830,]    4
 [831,]    8
 [832,]    7
 [833,]    4
 [834,]    5
 [835,]    5
 [836,]    5
 [837,]    5
 [838,]    3
 [839,]    5
 [840,]    3
 [841,]    6
 [842,]    8
 [843,]    9
 [844,]    4
 [845,]    8
 [846,]    2
 [847,]    6
 [848,]    4
 [849,]    4
 [850,]    3
 [851,]    7
 [852,]    9
 [853,]    8
 [854,]    4
 [855,]    7
 [856,]    4
 [857,]    8
 [858,]    4
 [859,]    5
 [860,]    8
 [861,]    7
 [862,]    5
 [863,]    6
 [864,]    7
 [865,]    9
 [866,]    8
 [867,]    7
 [868,]    8
 [869,]    6
 [870,]    6
 [871,]    7
 [872,]    8
 [873,]    9
 [874,]    4
 [875,]    0
 [876,]    1
 [877,]    5
 [878,]    6
 [879,]    5
 [880,]    2
 [881,]    1
 [882,]    9
 [883,]    9
 [884,]    7
 [885,]    5
 [886,]    5
 [887,]    5
 [888,]    5
 [889,]    6
 [890,]    5
 [891,]    7
 [892,]    8
 [893,]    8
 [894,]    6
 [895,]    5
 [896,]    5
 [897,]    7
 [898,]    8
 [899,]    4
 [900,]    5
 [901,]    6
 [902,]    7
 [903,]    4
 [904,]    4
 [905,]    7
 [906,]    9
 [907,]    7
 [908,]    6
 [909,]    6
 [910,]    8
 [911,]    8
 [912,]    7
 [913,]    5
 [914,]    8
 [915,]    5
 [916,]    8
 [917,]    7
 [918,]    7
 [919,]    8
 [920,]    4
 [921,]    7
 [922,]    3
 [923,]    2
 [924,]    4
 [925,]    6
 [926,]    8
 [927,]    4
 [928,]    3
 [929,]    5
 [930,]    5
 [931,]    7
 [932,]    4
 [933,]    6
 [934,]    6
 [935,]    3
 [936,]    6
 [937,]    4
 [938,]    7
 [939,]    7
 [940,]    9
 [941,]    5
 [942,]    5
 [943,]    5
 [944,]    9
 [945,]    4
 [946,]    4
 [947,]    1
 [948,]    4
 [949,]    6
 [950,]    1
 [951,]    6
 [952,]    6
 [953,]    9
 [954,]    5
 [955,]    3
 [956,]    3
 [957,]    6
 [958,]    8
 [959,]    7
 [960,]    5
 [961,]    6
 [962,]    7
 [963,]    5
 [964,]    6
 [965,]    4
 [966,]    4
 [967,]    4
 [968,]    4
 [969,]    4
 [970,]    6
 [971,]    3
 [972,]    8
 [973,]    6
 [974,]    2
 [975,]    4
 [976,]    5
 [977,]    7
 [978,]    8
 [979,]    8
 [980,]    9
 [981,]    7
 [982,]    8
 [983,]    6
 [984,]    7
 [985,]    7
 [986,]    4
 [987,]    4
 [988,]    5
 [989,]    5
 [990,]    2
 [991,]    5
 [992,]    5
 [993,]    6
 [994,]    7
 [995,]    9
 [996,]    8
 [997,]    7
 [998,]    3
 [999,]    4
[1000,]    2
[1001,]    7
[1002,]    8
[1003,]    8
[1004,]    5
[1005,]    7
[1006,]    9
[1007,]    5
[1008,]    9
[1009,]    8
[1010,]    7
[1011,]    8
[1012,]    4
[1013,]    6
[1014,]    9
[1015,]    6
[1016,]    5
[1017,]    5
[1018,]    6
[1019,]    5
[1020,]    6
[1021,]    1
[1022,]    3
[1023,]    6
[1024,]    5
[1025,]    6
[1026,]    3
[1027,]    7
[1028,]    8
[1029,]    7
[1030,]    8
[1031,]    9
[1032,]    7
[1033,]    6
[1034,]    7
[1035,]    6
[1036,]    3
[1037,]    6
[1038,]    8
[1039,]    5
[1040,]    7
[1041,]    4
[1042,]    6
[1043,]    6
[1044,]    5
[1045,]    5
[1046,]    3
[1047,]    7
[1048,]    3
[1049,]    9
[1050,]    8
[1051,]    6
[1052,]    8
[1053,]    5
[1054,]    7
[1055,]    8
[1056,]    6
[1057,]    3
[1058,]    3
[1059,]    3
[1060,]    4
[1061,]    7
[1062,]    8
[1063,]    9
[1064,]    8
[1065,]    7
[1066,]    8
[1067,]    7
[1068,]    9
[1069,]    7
[1070,]    7
[1071,]    4
[1072,]    6
[1073,]    8
[1074,]    5
[1075,]    4
[1076,]    7
[1077,]    8
[1078,]    6
[1079,]    2
[1080,]    6
[1081,]    7
[1082,]    7
[1083,]    5
[1084,]    5
[1085,]    4
[1086,]    8
[1087,]    6
[1088,]    6
[1089,]    4
[1090,]    4
[1091,]    6
[1092,]    6
[1093,]    8
[1094,]    7
[1095,]    7
[1096,]    7
[1097,]    8
[1098,]    6
[1099,]    2
[1100,]    5
[1101,]    4
[1102,]    7
[1103,]    6
[1104,]    8
[1105,]    9
[1106,]    5
[1107,]    4
[1108,]    4
[1109,]    5
[1110,]    6
[1111,]    7
[1112,]    6
[1113,]    7
[1114,]    5
[1115,]    7
[1116,]    7
[1117,]    8
[1118,]    7
[1119,]    4
[1120,]    5
[1121,]    7
[1122,]    5
[1123,]    6
[1124,]    6
[1125,]    8
[1126,]    6
[1127,]    4
[1128,]    7
[1129,]    2
[1130,]    3
[1131,]    5
[1132,]    3
[1133,]    6
[1134,]    7
[1135,]    5
[1136,]    4
[1137,]    7
[1138,]    5
[1139,]    4
[1140,]    3
[1141,]    3
[1142,]    3
[1143,]    7
[1144,]    9
[1145,]    1
[1146,]    7
[1147,]    9
[1148,]    8
[1149,]    9
[1150,]    8
[1151,]    6
[1152,]    4
[1153,]    7
[1154,]    7
[1155,]    9
[1156,]    6
[1157,]    6
[1158,]    4
[1159,]    6
[1160,]    6
[1161,]    6
[1162,]    5
[1163,]    5
[1164,]    6
[1165,]    7
[1166,]    4
[1167,]    5
[1168,]    4
[1169,]    6
[1170,]    7
[1171,]    5
[1172,]    6
[1173,]    3
[1174,]    2
[1175,]    3
[1176,]    9
[1177,]    9
[1178,]    7
[1179,]    6
[1180,]    8
[1181,]    4
[1182,]    4
[1183,]    4
[1184,]    7
[1185,]    3
[1186,]    7
[1187,]    5
[1188,]    8
[1189,]    4
[1190,]    9
[1191,]    7
[1192,]    4
[1193,]    8
[1194,]    3
[1195,]    6
[1196,]    4
[1197,]    7
[1198,]    8
[1199,]    6
[1200,]    8
[1201,]    6
[1202,]    5
[1203,]    9
[1204,]    2
[1205,]    3
[1206,]    5
[1207,]    5
[1208,]    7
[1209,]    2
[1210,]    6
[1211,]    7
[1212,]    9
[1213,]    7
[1214,]    8
[1215,]    8
[1216,]    8
[1217,]    7
[1218,]    8
[1219,]    7
[1220,]    7
[1221,]    3
[1222,]    8
[1223,]    4
[1224,]    9
[1225,]    4
[1226,]    5
[1227,]    7
[1228,]    6
[1229,]    7
[1230,]    6
[1231,]    4
[1232,]    5
[1233,]    4
[1234,]    4
[1235,]    2
[1236,]    4
[1237,]    4
[1238,]    4
[1239,]    6
[1240,]    5
[1241,]    7
[1242,]    7
[1243,]    7
[1244,]    9
[1245,]    5
[1246,]    8
[1247,]    9
[1248,]    8
[1249,]    8
[1250,]    7
[1251,]    7
[1252,]    8
[1253,]    9
[1254,]    6
[1255,]    7
[1256,]    8
[1257,]    6
[1258,]    4
[1259,]    6
[1260,]    6
[1261,]    7
[1262,]    6
[1263,]    6
[1264,]    6
[1265,]    4
[1266,]    3
[1267,]    3
[1268,]    4
[1269,]    3
[1270,]    5
[1271,]    8
[1272,]    6
[1273,]    5
[1274,]    7
[1275,]    7
[1276,]    8
[1277,]    5
[1278,]    7
[1279,]    7
[1280,]    3
[1281,]    9
[1282,]    5
[1283,]    2
[1284,]    6
[1285,]    4
[1286,]    3
[1287,]    6
[1288,]    4
[1289,]    4
[1290,]    6
[1291,]    4
[1292,]    8
[1293,]    6
[1294,]    6
[1295,]    7
[1296,]    9
[1297,]    9
[1298,]    7
[1299,]    6
[1300,]    2
[1301,]    7
[1302,]    6
[1303,]    7
[1304,]    7
[1305,]    7
[1306,]    9
[1307,]    7
[1308,]    5
[1309,]    7
[1310,]    7
[1311,]    5
[1312,]    7
[1313,]    7
[1314,]    7
[1315,]    6
[1316,]    5
[1317,]    3
[1318,]    5
[1319,]    6
[1320,]    6
[1321,]    8
[1322,]    6
[1323,]    6
[1324,]    5
[1325,]    5
[1326,]    5
[1327,]    5
[1328,]    8
[1329,]    8
[1330,]    5
[1331,]    6
[1332,]    6
[1333,]    4
[1334,]    5
[1335,]    7
[1336,]    7
[1337,]    4
[1338,]    6
[1339,]    6
[1340,]    5
[1341,]    7
[1342,]    8
[1343,]    7
[1344,]    8
[1345,]    8
[1346,]    4
[1347,]    4
[1348,]    3
[1349,]    6
[1350,]    3
[1351,]    4
[1352,]    6
[1353,]    6
[1354,]    7
[1355,]    9
[1356,]    8
[1357,]    4
[1358,]    6
[1359,]    6
[1360,]    4
[1361,]    6
[1362,]    5
[1363,]    3
[1364,]    3
[1365,]    6
[1366,]    3
[1367,]    9
[1368,]    7
[1369,]    6
[1370,]    3
[1371,]    6
[1372,]    3
[1373,]    6
[1374,]    5
[1375,]    5
[1376,]    7
[1377,]    7
[1378,]    5
[1379,]    2
[1380,]    5
[1381,]    4
[1382,]    4
[1383,]    3
[1384,]    6
[1385,]    6
[1386,]    8
[1387,]    5
[1388,]    5
[1389,]    6
[1390,]    6
[1391,]    7
[1392,]    7
[1393,]    4
[1394,]    3
[1395,]    7
[1396,]    6
[1397,]    6
[1398,]    7
[1399,]    6
[1400,]    0
[1401,]    3
[1402,]    3
[1403,]    5
[1404,]    5
[1405,]    5
[1406,]    7
[1407,]    9
[1408,]    6
[1409,]    1
[1410,]    9
[1411,]    9
[1412,]    5
[1413,]    5
[1414,]    5
[1415,]    7
[1416,]    7
[1417,]    6
[1418,]    9
[1419,]    7
[1420,]    7
[1421,]    7
[1422,]    4
[1423,]    8
[1424,]    5
[1425,]    6
[1426,]    7
[1427,]    3
[1428,]    7
[1429,]    3
[1430,]    5
[1431,]    3
[1432,]    5
[1433,]    3
[1434,]    6
[1435,]    1
[1436,]    5
[1437,]    2
[1438,]    5
[1439,]    3
[1440,]    5
[1441,]    1
[1442,]    4
[1443,]    5
[1444,]    6
[1445,]    1
[1446,]    5
[1447,]    8
[1448,]    8
[1449,]    5
[1450,]    5
[1451,]    3
[1452,]    8
[1453,]    9
[1454,]    6
[1455,]    7
[1456,]    7
[1457,]    9
[1458,]    7
[1459,]    7
[1460,]    6
[1461,]    6
[1462,]    7
[1463,]    7
[1464,]    3
[1465,]    6
[1466,]    5
[1467,]    7
[1468,]    8
[1469,]    6
[1470,]    6
[1471,]    7
[1472,]    8
[1473,]    7
[1474,]    8
[1475,]    2
[1476,]    6
[1477,]    6
[1478,]    4
[1479,]    8
[1480,]    5
[1481,]    5
[1482,]    2
[1483,]    8
[1484,]    9
[1485,]    9
[1486,]    8
[1487,]    3
[1488,]    7
[1489,]    6
[1490,]    6
[1491,]    6
[1492,]    8
[1493,]    8
[1494,]    5
[1495,]    8
[1496,]    7
[1497,]    7
[1498,]    7
[1499,]    4
[1500,]    5
[1501,]    2
[1502,]    8
[1503,]    6
[1504,]    5
[1505,]    8
[1506,]    8
[1507,]    7
[1508,]    7
[1509,]    6
[1510,]    4
[1511,]    5
[1512,]    6
[1513,]    6
[1514,]    6
[1515,]    6
[1516,]    7
[1517,]    3
[1518,]    5
[1519,]    4
[1520,]    3
[1521,]    4
[1522,]    3
[1523,]    8
[1524,]    8
[1525,]    9
[1526,]    6
[1527,]    7
[1528,]    4
[1529,]    6
[1530,]    7
[1531,]    4
[1532,]    4
[1533,]    3
[1534,]    8
[1535,]    9
[1536,]    6
[1537,]    7
[1538,]    8
[1539,]    3
[1540,]    6
[1541,]    6
[1542,]    6
[1543,]    7
[1544,]    5
[1545,]    5
[1546,]    5
[1547,]    5
[1548,]    5
[1549,]    6
[1550,]    8
[1551,]    9
[1552,]    7
[1553,]    3
[1554,]    7
[1555,]    7
[1556,]    5
[1557,]    4
[1558,]    6
[1559,]    7
[1560,]    2
[1561,]    8
[1562,]    6
[1563,]    4
[1564,]    4
[1565,]    5
[1566,]    3
[1567,]    6
[1568,]    3
[1569,]    3
[1570,]    7
[1571,]    1
[1572,]    5
[1573,]    9
[1574,]    9
[1575,]    9
[1576,]    8
[1577,]    3
[1578,]    3
[1579,]    3
[1580,]    2
[1581,]    5
[1582,]    4
[1583,]    5
[1584,]    7
[1585,]    7
[1586,]    8
[1587,]    3
[1588,]    4
[1589,]    7
[1590,]    6
[1591,]    6
[1592,]    4
[1593,]    3
[1594,]    8
[1595,]    7
[1596,]    5
[1597,]    5
[1598,]    7
[1599,]    8
[1600,]    6
[1601,]    7
[1602,]    5
[1603,]    9
[1604,]    5
[1605,]    7
[1606,]    8
[1607,]    8
[1608,]    7
[1609,]    7
[1610,]    7
[1611,]    5
[1612,]    5
[1613,]    7
[1614,]    4
[1615,]    4
[1616,]    3
[1617,]    5
[1618,]    7
[1619,]    6
[1620,]    5
[1621,]    5
[1622,]    9
[1623,]    5
[1624,]    7
[1625,]    7
[1626,]    4
[1627,]    5
[1628,]    6
[1629,]    7
[1630,]    4
[1631,]    7
[1632,]    6
[1633,]    4
[1634,]    9
[1635,]    5
[1636,]    7
[1637,]    5
[1638,]    7
[1639,]    2
[1640,]    6
[1641,]    3
[1642,]    6
[1643,]    5
[1644,]    6
[1645,]    5
[1646,]    7
[1647,]    8
[1648,]    7
[1649,]    7
[1650,]    8
[1651,]    5
[1652,]    6
[1653,]    6
[1654,]    6
[1655,]    8
[1656,]    7
[1657,]    4
[1658,]    6
[1659,]    5
[1660,]    7
[1661,]    5
[1662,]    7
[1663,]    4
[1664,]    5
[1665,]    3
[1666,]    5
[1667,]    4
[1668,]    7
[1669,]    8
[1670,]    6
[1671,]    8
[1672,]    8
[1673,]    7
[1674,]    6
[1675,]    5
[1676,]    4
[1677,]    7
[1678,]    4
[1679,]    5
[1680,]    2
[1681,]    6
[1682,]    6
[1683,]    5
[1684,]    6
[1685,]    3
[1686,]    6
[1687,]    6
[1688,]    8
[1689,]    6
[1690,]    4
[1691,]    4
[1692,]    7
[1693,]    4
[1694,]    5
[1695,]    8
[1696,]    4
[1697,]    5
[1698,]    4
[1699,]    7
[1700,]    5
[1701,]    9
[1702,]    4
[1703,]    7
[1704,]    5
[1705,]    7
[1706,]    3
[1707,]    8
[1708,]    5
[1709,]    8
[1710,]    5
[1711,]    5
[1712,]    4
[1713,]    9
[1714,]    6
[1715,]    4
[1716,]    7
[1717,]    3
[1718,]    5
[1719,]    6
[1720,]    7
[1721,]    5
[1722,]    8
[1723,]    4
[1724,]    4
[1725,]    4
[1726,]    6
[1727,]    3
[1728,]    3
[1729,]    5
[1730,]    3
[1731,]    0
[1732,]    2
[1733,]    3
[1734,]    4
[1735,]    5
[1736,]    3
[1737,]    5
[1738,]    5
[1739,]    5
[1740,]    4
[1741,]    8
[1742,]    9
[1743,]    8
[1744,]    5
[1745,]    8
[1746,]    7
[1747,]    7
[1748,]    8
[1749,]    6
[1750,]    4
[1751,]    6
[1752,]    5
[1753,]    7
[1754,]    7
[1755,]    6
[1756,]    5
[1757,]    5
[1758,]    2
[1759,]    1
[1760,]    3
[1761,]    4
[1762,]    9
[1763,]    7
[1764,]    5
[1765,]    8
[1766,]    6
[1767,]    9
[1768,]    6
[1769,]    5
[1770,]    6
[1771,]    6
[1772,]    3
[1773,]    3
[1774,]    4
[1775,]    5
[1776,]    5
[1777,]    4
[1778,]    8
[1779,]    7
[1780,]    6
[1781,]    8
[1782,]    6
[1783,]    7
[1784,]    8
[1785,]    8
[1786,]    5
[1787,]    5
[1788,]    5
[1789,]    4
[1790,]    6
[1791,]    6
[1792,]    8
[1793,]    9
[1794,]    6
[1795,]    6
[1796,]    4
[1797,]    6
[1798,]    4
[1799,]    2
[1800,]    4
[1801,]    4
[1802,]    6
[1803,]    5
[1804,]    6
[1805,]    6
[1806,]    9
[1807,]    4
[1808,]    4
[1809,]    4
[1810,]    8
[1811,]    5
[1812,]    6
[1813,]    7
[1814,]    6
[1815,]    6
[1816,]    6
[1817,]    6
[1818,]    4
[1819,]    4
[1820,]    3
[1821,]    4
[1822,]    4
[1823,]    2
[1824,]    5
[1825,]    4
[1826,]    7
[1827,]    5
[1828,]    4
[1829,]    8
[1830,]    9
[1831,]    9
[1832,]    9
[1833,]    8
[1834,]    4
[1835,]    9
[1836,]    2
[1837,]    8
[1838,]    5
[1839,]    5
[1840,]    8
[1841,]    7
[1842,]    3
[1843,]    3
[1844,]    5
[1845,]    4
[1846,]    6
[1847,]    8
[1848,]    6
[1849,]    7
[1850,]    5
[1851,]    7
[1852,]    5
[1853,]    5
[1854,]    7
[1855,]    6
[1856,]    8
[1857,]    3
[1858,]    7
[1859,]    5
[1860,]    4
[1861,]    4
[1862,]    2
[1863,]    5
[1864,]    4
[1865,]    4
[1866,]    6
[1867,]    3
[1868,]    4
[1869,]    7
[1870,]    9
[1871,]    7
[1872,]    3
[1873,]    6
[1874,]    5
[1875,]    6
[1876,]    6
[1877,]    7
[1878,]    2
[1879,]    3
[1880,]    5
[1881,]    6
[1882,]    9
[1883,]    6
[1884,]    7
[1885,]    6
[1886,]    7
[1887,]    3
[1888,]    3
[1889,]    7
[1890,]    7
[1891,]    5
[1892,]    4
[1893,]    7
[1894,]    5
[1895,]    9
[1896,]    6
[1897,]    3
[1898,]    8
[1899,]    4
[1900,]    5
[1901,]    4
[1902,]    5
[1903,]    5
[1904,]    7
[1905,]    5
[1906,]    8
[1907,]    5
[1908,]    7
[1909,]    9
[1910,]    8
[1911,]    6
[1912,]    5
[1913,]    5
[1914,]    5
[1915,]    5
[1916,]    7
[1917,]    7
[1918,]    8
[1919,]    4
[1920,]    7
[1921,]    3
[1922,]    5
[1923,]    6
[1924,]    6
[1925,]    9
[1926,]    8
[1927,]    6
[1928,]    3
[1929,]    9
[1930,]    8
[1931,]    7
[1932,]    4
[1933,]    5
[1934,]    5
[1935,]    5
[1936,]    6
[1937,]    4
[1938,]    0
[1939,]    7
[1940,]    0
[1941,]    4
[1942,]    7
[1943,]    4
[1944,]    1
[1945,]    9
[1946,]    5
[1947,]    5
[1948,]    4
[1949,]    4
[1950,]    1
[1951,]    8
[1952,]    6
[1953,]    6
[1954,]    7
[1955,]    5
[1956,]    7
[1957,]    6
[1958,]    8
[1959,]    9
[1960,]    6
[1961,]    8
[1962,]    3
[1963,]    6
[1964,]    3
[1965,]    3
[1966,]    4
[1967,]    4
[1968,]    4
[1969,]    5
[1970,]    5
[1971,]    5
[1972,]    5
[1973,]    5
[1974,]    7
[1975,]    5
[1976,]    8
[1977,]    6
[1978,]    6
[1979,]    3
[1980,]    6
[1981,]    5
[1982,]    4
[1983,]    5
[1984,]    7
[1985,]    6
[1986,]    7
[1987,]    6
[1988,]    8
[1989,]    6
[1990,]    3
[1991,]    7
[1992,]    4
[1993,]    4
[1994,]    1
[1995,]    3
[1996,]    4
[1997,]    1
[1998,]    4
[1999,]    3
[2000,]    7
[2001,]    7
[2002,]    3
[2003,]    1
[2004,]    4
[2005,]    7
[2006,]    6
[2007,]    8
[2008,]    1
[2009,]    2
[2010,]    5
[2011,]    7
[2012,]    3
[2013,]    8
[2014,]    7
[2015,]    7
[2016,]    4
[2017,]    5
[2018,]    1
[2019,]    7
[2020,]    4
[2021,]    7
[2022,]    6
[2023,]    4
[2024,]    6
[2025,]    6
[2026,]    9
[2027,]    7
[2028,]    4
[2029,]    7
[2030,]    7
[2031,]    5
[2032,]    8
[2033,]    6
[2034,]    5
[2035,]    6
[2036,]    4
[2037,]    7
[2038,]    6
[2039,]    5
[2040,]    5
[2041,]    5
[2042,]    5
[2043,]    4
[2044,]    8
[2045,]    4
[2046,]    8
[2047,]    7
[2048,]    6
[2049,]    3
[2050,]    6
[2051,]    5
[2052,]    6
[2053,]    5
[2054,]    5
[2055,]    5
[2056,]    5
[2057,]    5
[2058,]    4
[2059,]    6
[2060,]    1
[2061,]    8
[2062,]    5
[2063,]    3
[2064,]    8
[2065,]    9
[2066,]    7
[2067,]    7
[2068,]    7
[2069,]    7
[2070,]    7
[2071,]    5
[2072,]    3
[2073,]    8
[2074,]    5
[2075,]    6
[2076,]    3
[2077,]    2
[2078,]    3
[2079,]    5
[2080,]    3
[2081,]    1
[2082,]    6
[2083,]    8
[2084,]    7
[2085,]    5
[2086,]    6
[2087,]    2
[2088,]    1
[2089,]    6
[2090,]    5
[2091,]    6
[2092,]    7
[2093,]    5
[2094,]    2
[2095,]    2
[2096,]    3
[2097,]    5
[2098,]    8
[2099,]    6
[2100,]    7
[2101,]    7
[2102,]    4
[2103,]    5
[2104,]    3
[2105,]    6
[2106,]    9
[2107,]    8
[2108,]    7
[2109,]    8
[2110,]    5
[2111,]    6
[2112,]    9
[2113,]    9
[2114,]    7
[2115,]    9
[2116,]    7
[2117,]    7
[2118,]    6
[2119,]    4
[2120,]    8
[2121,]    6
[2122,]    2
[2123,]    3
[2124,]    4
[2125,]    6
[2126,]    6
[2127,]    8
[2128,]    4
[2129,]    2
[2130,]    2
[2131,]    4
[2132,]    6
[2133,]    6
[2134,]    6
[2135,]    8
[2136,]    4
[2137,]    7
[2138,]    4
[2139,]    8
[2140,]    4
[2141,]    7
[2142,]    7
[2143,]    5
[2144,]    6
[2145,]    7
[2146,]    6
[2147,]    3
[2148,]    6
[2149,]    2
[2150,]    5
[2151,]    4
[2152,]    9
[2153,]    9
[2154,]    7
[2155,]    6
[2156,]    3
[2157,]    4
[2158,]    8
[2159,]    7
[2160,]    6
[2161,]    6
[2162,]    5
[2163,]    6
[2164,]    7
[2165,]    7
[2166,]    7
[2167,]    4
[2168,]    1
[2169,]    5
[2170,]    6
[2171,]    7
[2172,]    6
[2173,]    4
[2174,]    6
[2175,]    7
[2176,]    4
[2177,]    5
[2178,]    8
[2179,]    6
[2180,]    6
[2181,]    3
[2182,]    4
[2183,]    5
[2184,]    5
[2185,]    6
[2186,]    4
[2187,]    7
[2188,]    7
[2189,]    8
[2190,]    9
[2191,]    6
[2192,]    7
[2193,]    6
[2194,]    5
[2195,]    7
[2196,]    8
[2197,]    6
[2198,]    7
[2199,]    8
[2200,]    8
[2201,]    3
[2202,]    5
[2203,]    7
[2204,]    8
[2205,]    7
[2206,]    5
[2207,]    8
[2208,]    7
[2209,]    5
[2210,]    5
[2211,]    7
[2212,]    6
[2213,]    5
[2214,]    5
[2215,]    3
[2216,]    2
[2217,]    4
[2218,]    4
[2219,]    5
[2220,]    6
[2221,]    5
[2222,]    5
[2223,]    3
[2224,]    3
[2225,]    9
[2226,]    5
[2227,]    5
[2228,]    4
[2229,]    5
[2230,]    6
[2231,]    6
[2232,]    7
[2233,]    6
[2234,]    6
[2235,]    5
[2236,]    6
[2237,]    9
[2238,]    6
[2239,]    5
[2240,]    7
[2241,]    9
[2242,]    6
[2243,]    5
[2244,]    7
[2245,]    8
[2246,]    7
[2247,]    3
[2248,]    8
[2249,]    7
[2250,]    8
[2251,]    7
[2252,]    7
[2253,]    4
[2254,]    7
[2255,]    7
[2256,]    5
[2257,]    3
[2258,]    7
[2259,]    5
[2260,]    8
[2261,]    5
[2262,]    5
[2263,]    7
[2264,]    7
[2265,]    6
[2266,]    6
[2267,]    4
[2268,]    5
[2269,]    9
[2270,]    7
[2271,]    7
[2272,]    6
[2273,]    4
[2274,]    7
[2275,]    6
[2276,]    6
[2277,]    8
[2278,]    7
[2279,]    5
[2280,]    5
[2281,]    7
[2282,]    7
[2283,]    5
[2284,]    8
[2285,]    4
[2286,]    6
[2287,]    7
[2288,]    6
[2289,]    8
[2290,]    3
[2291,]    6
[2292,]    9
[2293,]    5
[2294,]    7
[2295,]    9
[2296,]    9
[2297,]    5
[2298,]    6
[2299,]    5
[2300,]    4
[2301,]    6
[2302,]    7
[2303,]    7
[2304,]    3
[2305,]    3
[2306,]    8
[2307,]    9
[2308,]    3
[2309,]    7
[2310,]    8
[2311,]    7
[2312,]    4
[2313,]    3
[2314,]    5
[2315,]    2
[2316,]    5
[2317,]    3
[2318,]    7
[2319,]    4
[2320,]    4
[2321,]    4
[2322,]    7
[2323,]    5
[2324,]    6
[2325,]    7
[2326,]    8
[2327,]    5
[2328,]    6
[2329,]    7
[2330,]    4
[2331,]    6
[2332,]    7
[2333,]    5
[2334,]    7
[2335,]    6
[2336,]    3
[2337,]    5
[2338,]    6
[2339,]    4
[2340,]    5
[2341,]    8
[2342,]    8
[2343,]    8
[2344,]    9
[2345,]    7
[2346,]    4
[2347,]    8
[2348,]    8
[2349,]    6
[2350,]    5
[2351,]    6
[2352,]    9
[2353,]    3
[2354,]    6
[2355,]    4
[2356,]    2
[2357,]    6
[2358,]    7
[2359,]    3
[2360,]    4
[2361,]    5
[2362,]    7
[2363,]    7
[2364,]    8
[2365,]    5
[2366,]    2
[2367,]    3
[2368,]    7
[2369,]    2
[2370,]    9
[2371,]    7
[2372,]    6
[2373,]    7
[2374,]    5
[2375,]    6
[2376,]    5
[2377,]    5
[2378,]    8
[2379,]    6
[2380,]    3
[2381,]    6
[2382,]    6
[2383,]    8
[2384,]    8
[2385,]    6
[2386,]    4
[2387,]    3
[2388,]    6
[2389,]    0
[2390,]    2
[2391,]    6
[2392,]    6
[2393,]    3
[2394,]    7
[2395,]    4
[2396,]    7
[2397,]    4
[2398,]    7
[2399,]    7
[2400,]    4
[2401,]    8
[2402,]    8
[2403,]    5
[2404,]    4
[2405,]    8
[2406,]    8
[2407,]    7
[2408,]    6
[2409,]    7
[2410,]    5
[2411,]    6
[2412,]    4
[2413,]    5
[2414,]    9
[2415,]    4
[2416,]    6
[2417,]    6
[2418,]    8
[2419,]    7
[2420,]    7
[2421,]    3
[2422,]    4
[2423,]    7
[2424,]    4
[2425,]    5
[2426,]    4
[2427,]    7
[2428,]    7
[2429,]    4
[2430,]    7
[2431,]    4
[2432,]    4
[2433,]    4
[2434,]    5
[2435,]    3
[2436,]    7
[2437,]    7
[2438,]    6
[2439,]    6
[2440,]    7
[2441,]    5
[2442,]    3
[2443,]    4
[2444,]    4
[2445,]    7
[2446,]    3
[2447,]    4
[2448,]    4
[2449,]    7
[2450,]    6
[2451,]    4
[2452,]    5
[2453,]    4
[2454,]    7
[2455,]    7
[2456,]    6
[2457,]    7
[2458,]    2
[2459,]    8
[2460,]    6
[2461,]    6
[2462,]    8
[2463,]    6
[2464,]    8
[2465,]    1
[2466,]    4
[2467,]    8
[2468,]    6
[2469,]    7
[2470,]    4
[2471,]    8
[2472,]    6
[2473,]    5
[2474,]    5
[2475,]    6
[2476,]    8
[2477,]    6
[2478,]    9
[2479,]    9
[2480,]    7
[2481,]    9
[2482,]    5
[2483,]    6
[2484,]    3
[2485,]    5
[2486,]    7
[2487,]    7
[2488,]    7
[2489,]    5
[2490,]    7
[2491,]    6
[2492,]    5
[2493,]    4
[2494,]    8
[2495,]    8
[2496,]    7
[2497,]    4
[2498,]    6
[2499,]    8
[2500,]    5
[2501,]    7
[2502,]    8
[2503,]    7
[2504,]    7
[2505,]    7
[2506,]    6
[2507,]    7
[2508,]    7
[2509,]    6
[2510,]    9
[2511,]    7
[2512,]    5
[2513,]    3
[2514,]    8
[2515,]    6
[2516,]    9
[2517,]    3
[2518,]    9
[2519,]    9
[2520,]    7
[2521,]    6
[2522,]    7
[2523,]    8
[2524,]    7
[2525,]    6
[2526,]    4
[2527,]    5
[2528,]    7
[2529,]    6
[2530,]    7
[2531,]    9
[2532,]    8
[2533,]    5
[2534,]    7
[2535,]    5
[2536,]    8
[2537,]    5
[2538,]    7
[2539,]    4
[2540,]    6
[2541,]    3
[2542,]    3
[2543,]    3
[2544,]    2
[2545,]    4
[2546,]    6
[2547,]    5
[2548,]    6
[2549,]    2
[2550,]    7
[2551,]    6
[2552,]    6
[2553,]    8
[2554,]    2
[2555,]    6
[2556,]    6
[2557,]    8
[2558,]    6
[2559,]    6
[2560,]    6
[2561,]    9
[2562,]    7
[2563,]    5
[2564,]    8
[2565,]    8
[2566,]    6
[2567,]    5
[2568,]    7
[2569,]    1
[2570,]    6
[2571,]    8
[2572,]    4
[2573,]    7
[2574,]    6
[2575,]    5
[2576,]    5
[2577,]    4
[2578,]    6
[2579,]    6
[2580,]    6
[2581,]    6
[2582,]    6
[2583,]    8
[2584,]    9
[2585,]    5
[2586,]    8
[2587,]    8
[2588,]    8
[2589,]    6
[2590,]    4
[2591,]    7
[2592,]    7
[2593,]    6
[2594,]    4
[2595,]    4
[2596,]    5
[2597,]    8
[2598,]    7
[2599,]    8
[2600,]    1
[2601,]    4
[2602,]    5
[2603,]    8
[2604,]    8
[2605,]    5
[2606,]    5
[2607,]    5
[2608,]    5
[2609,]    7
[2610,]    7
[2611,]    6
[2612,]    6
[2613,]    4
[2614,]    8
[2615,]    3
[2616,]    5
[2617,]    5
[2618,]    6
[2619,]    6
[2620,]    7
[2621,]    6
[2622,]    4
[2623,]    3
[2624,]    4
[2625,]    4
[2626,]    5
[2627,]    3
[2628,]    6
[2629,]    8
[2630,]    7
[2631,]    8
[2632,]    6
[2633,]    3
[2634,]    3
[2635,]    4
[2636,]    8
[2637,]    7
[2638,]    7
[2639,]    8
[2640,]    7
[2641,]    8
[2642,]    7
[2643,]    6
[2644,]    8
[2645,]    7
[2646,]    7
[2647,]    6
[2648,]    4
[2649,]    8
[2650,]    4
[2651,]    7
[2652,]    3
[2653,]    5
[2654,]    1
[2655,]    1
[2656,]    6
[2657,]    4
[2658,]    5
[2659,]    3
[2660,]    6
[2661,]    6
[2662,]    4
[2663,]    9
[2664,]    8
[2665,]    6
[2666,]    7
[2667,]    4
[2668,]    6
[2669,]    7
[2670,]    8
[2671,]    7
[2672,]    8
[2673,]    2
[2674,]    6
[2675,]    7
[2676,]    5
[2677,]    4
[2678,]    4
[2679,]    7
[2680,]    6
[2681,]    9
[2682,]    7
[2683,]    9
[2684,]    5
[2685,]    6
[2686,]    6
[2687,]    9
[2688,]    6
[2689,]    6
[2690,]    6
[2691,]    9
[2692,]    9
[2693,]    5
[2694,]    6
[2695,]    5
[2696,]    4
[2697,]    2
[2698,]    3
[2699,]    5
[2700,]    3
[2701,]    7
[2702,]    2
[2703,]    4
[2704,]    3
[2705,]    4
[2706,]    0
[2707,]    5
[2708,]    6
[2709,]    7
[2710,]    8
[2711,]    1
[2712,]    5
[2713,]    6
[2714,]    5
[2715,]    7
[2716,]    3
[2717,]    7
[2718,]    8
[2719,]    4
[2720,]    7
[2721,]    6
[2722,]    7
[2723,]    5
[2724,]    2
[2725,]    7
[2726,]    8
[2727,]    6
[2728,]    7
[2729,]    7
[2730,]    6
[2731,]    5
[2732,]    7
[2733,]    6
[2734,]    4
[2735,]    8
[2736,]    6
[2737,]    4
[2738,]    5
[2739,]    7
[2740,]    5
[2741,]    6
[2742,]    6
[2743,]    4
[2744,]    4
[2745,]    1
[2746,]    8
[2747,]    9
[2748,]    8
[2749,]    9
[2750,]    8
[2751,]    8
[2752,]    3
[2753,]    5
[2754,]    6
[2755,]    7
[2756,]    9
[2757,]    7
[2758,]    6
[2759,]    9
[2760,]    7
[2761,]    6
[2762,]    7
[2763,]    9
[2764,]    9
[2765,]    8
[2766,]    8
[2767,]    4
[2768,]    3
[2769,]    3
[2770,]    5
[2771,]    5
[2772,]    7
[2773,]    4
[2774,]    8
[2775,]    4
[2776,]    9
[2777,]    8
[2778,]    7
[2779,]    8
[2780,]    4
[2781,]    8
[2782,]    5
[2783,]    7
[2784,]    3
[2785,]    2
[2786,]    7
[2787,]    4
[2788,]    6
[2789,]    3
[2790,]    8
[2791,]    7
[2792,]    6
[2793,]    2
[2794,]    0
[2795,]    3
[2796,]    4
[2797,]    8
[2798,]    6
[2799,]    7
[2800,]    8
[2801,]    7
[2802,]    5
[2803,]    5
[2804,]    7
[2805,]    8
[2806,]    4
[2807,]    7
[2808,]    5
[2809,]    4
[2810,]    5
[2811,]    4
[2812,]    5
[2813,]    7
[2814,]    7
[2815,]    9
[2816,]    7
[2817,]    4
[2818,]    8
[2819,]    6
[2820,]    7
[2821,]    6
[2822,]    6
[2823,]    8
[2824,]    6
[2825,]    7
[2826,]    6
[2827,]    2
[2828,]    6
[2829,]    1
[2830,]    4
[2831,]    5
[2832,]    9
[2833,]    9
[2834,]    8
[2835,]    9
[2836,]    6
[2837,]    6
[2838,]    7
[2839,]    6
[2840,]    8
[2841,]    8
[2842,]    8
[2843,]    4
[2844,]    3
[2845,]    3
[2846,]    4
[2847,]    5
[2848,]    6
[2849,]    6
[2850,]    5
[2851,]    7
[2852,]    5
[2853,]    8
[2854,]    5
[2855,]    5
[2856,]    3
[2857,]    3
[2858,]    3
[2859,]    2
[2860,]    7
[2861,]    6
[2862,]    5
[2863,]    3
[2864,]    6
[2865,]    3
[2866,]    5
[2867,]    4
[2868,]    7
[2869,]    6
[2870,]    5
[2871,]    4
[2872,]    6
[2873,]    7
[2874,]    3
[2875,]    4
[2876,]    5
[2877,]    3
[2878,]    4
[2879,]    7
[2880,]    7
[2881,]    4
[2882,]    4
[2883,]    8
[2884,]    8
[2885,]    9
[2886,]    6
[2887,]    8
[2888,]    6
[2889,]    4
[2890,]    6
[2891,]    7
[2892,]    4
[2893,]    2
[2894,]    6
[2895,]    9
[2896,]    8
[2897,]    4
[2898,]    5
[2899,]    4
[2900,]    6
[2901,]    7
[2902,]    6
[2903,]    4
[2904,]    4
[2905,]    8
[2906,]    7
[2907,]    8
[2908,]    3
[2909,]    2
[2910,]    3
[2911,]    2
[2912,]    2
[2913,]    2
[2914,]    1
[2915,]    5
[2916,]    8
[2917,]    7
[2918,]    5
[2919,]    8
[2920,]    6
[2921,]    7
[2922,]    8
[2923,]    9
[2924,]    5
[2925,]    7
[2926,]    7
[2927,]    3
[2928,]    2
[2929,]    3
[2930,]    6
[2931,]    8
[2932,]    7
[2933,]    8
[2934,]    7
[2935,]    7
[2936,]    6
[2937,]    6
[2938,]    5
[2939,]    4
[2940,]    8
[2941,]    8
[2942,]    5
[2943,]    6
[2944,]    7
[2945,]    6
[2946,]    7
[2947,]    5
[2948,]    7
[2949,]    6
[2950,]    7
[2951,]    3
[2952,]    4
[2953,]    8
[2954,]    6
[2955,]    7
[2956,]    5
[2957,]    5
[2958,]    6
[2959,]    4
[2960,]    6
[2961,]    5
[2962,]    5
[2963,]    5
[2964,]    9
[2965,]    8
[2966,]    8
[2967,]    9
[2968,]    8
[2969,]    8
[2970,]    8
[2971,]    6
[2972,]    6
[2973,]    5
[2974,]    8
[2975,]    8
[2976,]    8
[2977,]    2
[2978,]    6
[2979,]    7
[2980,]    3
[2981,]    7
[2982,]    2
[2983,]    7
[2984,]    4
[2985,]    7
[2986,]    7
[2987,]    8
[2988,]    6
[2989,]    6
[2990,]    6
[2991,]    3
[2992,]    5
[2993,]    7
[2994,]    5
[2995,]    4
[2996,]    3
[2997,]    5
[2998,]    7
[2999,]    5
[3000,]    4
[3001,]    6
[3002,]    6
[3003,]    7
[3004,]    8
[3005,]    8
[3006,]    4
[3007,]    6
[3008,]    7
[3009,]    3
[3010,]    6
[3011,]    3
[3012,]    6
[3013,]    8
[3014,]    8
[3015,]    7
[3016,]    6
[3017,]    5
[3018,]    7
[3019,]    6
[3020,]    6
[3021,]    4
[3022,]    7
[3023,]    7
[3024,]    3
[3025,]    4
[3026,]    4
[3027,]    4
[3028,]    6
[3029,]    8
[3030,]    9
[3031,]    8
[3032,]    8
[3033,]    7
[3034,]    6
[3035,]    2
[3036,]    4
[3037,]    8
[3038,]    6
[3039,]    5
[3040,]    6
[3041,]    5
[3042,]    3
[3043,]    5
[3044,]    4
[3045,]    5
[3046,]    3
[3047,]    7
[3048,]    7
[3049,]    7
[3050,]    7
[3051,]    5
[3052,]    7
[3053,]    6
[3054,]    3
[3055,]    7
[3056,]    4
[3057,]    4
[3058,]    7
[3059,]    4
[3060,]    2
[3061,]    3
[3062,]    4
[3063,]    8
[3064,]    3
[3065,]    3
[3066,]    3
[3067,]    6
[3068,]    4
[3069,]    9
[3070,]    9
[3071,]    7
[3072,]    5
[3073,]    8
[3074,]    5
[3075,]    6
[3076,]    7
[3077,]    4
[3078,]    6
[3079,]    6
[3080,]    5
[3081,]    7
[3082,]    4
[3083,]    4
[3084,]    4
[3085,]    3
[3086,]    8
[3087,]    3
[3088,]    7
[3089,]    7
[3090,]    9
[3091,]    6
[3092,]    6
[3093,]    7
[3094,]    3
[3095,]    5
[3096,]    6
[3097,]    6
[3098,]    7
[3099,]    6
[3100,]    5
[3101,]    5
[3102,]    1
[3103,]    8
[3104,]    7
[3105,]    8
[3106,]    4
[3107,]    7
[3108,]    7
[3109,]    8
[3110,]    5
[3111,]    7
[3112,]    4
[3113,]    6
[3114,]    6
[3115,]    6
[3116,]    6
[3117,]    4
[3118,]    7
[3119,]    3
[3120,]    5
[3121,]    5
[3122,]    8
[3123,]    4
[3124,]    2
[3125,]    1
[3126,]    7
[3127,]    7
[3128,]    7
[3129,]    0
[3130,]    5
[3131,]    4
[3132,]    3
[3133,]    3
[3134,]    4
[3135,]    3
[3136,]    7
[3137,]    5
[3138,]    5
[3139,]    6
[3140,]    7
[3141,]    4
[3142,]    5
[3143,]    3
[3144,]    5
[3145,]    3
[3146,]    2
[3147,]    3
[3148,]    2
[3149,]    9
[3150,]    4
[3151,]    9
[3152,]    6
[3153,]    6
[3154,]    4
[3155,]    4
[3156,]    8
[3157,]    2
[3158,]    4
[3159,]    3
[3160,]    5
[3161,]    4
[3162,]    6
[3163,]    6
[3164,]    2
[3165,]    5
[3166,]    8
[3167,]    3
[3168,]    5
[3169,]    6
[3170,]    6
[3171,]    5
[3172,]    6
[3173,]    9
[3174,]    7
[3175,]    4
[3176,]    7
[3177,]    8
[3178,]    3
[3179,]    5
[3180,]    9
[3181,]    4
[3182,]    7
[3183,]    8
[3184,]    7
[3185,]    6
[3186,]    4
[3187,]    7
[3188,]    6
[3189,]    6
[3190,]    5
[3191,]    3
[3192,]    4
[3193,]    5
[3194,]    4
[3195,]    1
[3196,]    4
[3197,]    2
[3198,]    5
[3199,]    7
[3200,]    6
[3201,]    3
[3202,]    6
[3203,]    8
[3204,]    6
[3205,]    7
[3206,]    5
[3207,]    3
[3208,]    5
[3209,]    5
[3210,]    1
[3211,]    4
[3212,]    4
[3213,]    5
[3214,]    5
[3215,]    7
[3216,]    5
[3217,]    8
[3218,]    6
[3219,]    7
[3220,]    3
[3221,]    9
[3222,]    5
[3223,]    4
[3224,]    4
[3225,]    7
[3226,]    6
[3227,]    7
[3228,]    6
[3229,]    5
[3230,]    5
[3231,]    2
[3232,]    3
[3233,]    9
[3234,]    8
[3235,]    7
[3236,]    6
[3237,]    4
[3238,]    7
[3239,]    5
[3240,]    6
[3241,]    5
[3242,]    5
[3243,]    5
[3244,]    7
[3245,]    5
[3246,]    5
[3247,]    5
[3248,]    6
[3249,]    5
[3250,]    8
[3251,]    2
[3252,]    8
[3253,]    5
[3254,]    5
[3255,]    8
[3256,]    4
[3257,]    8
[3258,]    5
[3259,]    5
[3260,]    7
[3261,]    6
[3262,]    6
[3263,]    5
[3264,]    8
[3265,]    7
[3266,]    3
[3267,]    6
[3268,]    7
[3269,]    4
[3270,]    5
[3271,]    3
[3272,]    3
[3273,]    5
[3274,]    7
[3275,]    8
[3276,]    7
[3277,]    7
[3278,]    3
[3279,]    3
[3280,]    6
[3281,]    3
[3282,]    7
[3283,]    5
[3284,]    9
[3285,]    7
[3286,]    4
[3287,]    5
[3288,]    4
[3289,]    7
[3290,]    4
[3291,]    9
[3292,]    8
[3293,]    6
[3294,]    6
[3295,]    6
[3296,]    3
[3297,]    5
[3298,]    7
[3299,]    6
[3300,]    6
[3301,]    8
[3302,]    5
[3303,]    6
[3304,]    6
[3305,]    6
[3306,]    5
[3307,]    5
[3308,]    9
[3309,]    8
[3310,]    5
[3311,]    6
[3312,]    2
[3313,]    3
[3314,]    5
[3315,]    8
[3316,]    5
[3317,]    5
[3318,]    5
[3319,]    8
[3320,]    7
[3321,]    5
[3322,]    2
[3323,]    7
[3324,]    4
[3325,]    5
[3326,]    4
[3327,]    5
[3328,]    6
[3329,]    6
[3330,]    6
[3331,]    7
[3332,]    4
[3333,]    7
[3334,]    5
[3335,]    6
[3336,]    5
[3337,]    7
[3338,]    9
[3339,]    9
[3340,]    7
[3341,]    7
[3342,]    3
[3343,]    6
[3344,]    6
[3345,]    3
[3346,]    3
[3347,]    5
[3348,]    4
[3349,]    7
[3350,]    7
[3351,]    6
[3352,]    5
[3353,]    3
[3354,]    3
[3355,]    4
[3356,]    4
[3357,]    4
[3358,]    6
[3359,]    7
[3360,]    3
[3361,]    7
[3362,]    5
[3363,]    5
[3364,]    5
[3365,]    4
[3366,]    3
[3367,]    9
[3368,]    9
[3369,]    8
[3370,]    8
[3371,]    9
[3372,]    7
[3373,]    8
[3374,]    8
[3375,]    7
[3376,]    9
[3377,]    6
[3378,]    6
[3379,]    4
[3380,]    6
[3381,]    9
[3382,]    7
[3383,]    5
[3384,]    7
[3385,]    5
[3386,]    4
[3387,]    8
[3388,]    8
[3389,]    5
[3390,]    6
[3391,]    7
[3392,]    5
[3393,]    4
[3394,]    9
[3395,]    9
[3396,]    7
[3397,]    8
[3398,]    9
[3399,]    9
[3400,]    5
[3401,]    6
[3402,]    6
[3403,]    9
[3404,]    8
[3405,]    5
[3406,]    5
[3407,]    8
[3408,]    4
[3409,]    6
[3410,]    6
[3411,]    8
[3412,]    5
[3413,]    6
[3414,]    6
[3415,]    0
[3416,]    4
[3417,]    7
[3418,]    6
[3419,]    6
[3420,]    9
[3421,]    8
[3422,]    6
[3423,]    6
[3424,]    9
[3425,]    7
[3426,]    6
[3427,]    7
[3428,]    8
[3429,]    6
[3430,]    8
[3431,]    8
[3432,]    7
[3433,]    9
[3434,]    7
[3435,]    8
[3436,]    9
[3437,]    3
[3438,]    7
[3439,]    6
[3440,]    4
[3441,]    5
[3442,]    4
[3443,]    7
[3444,]    8
[3445,]    7
[3446,]    8
[3447,]    7
[3448,]    5
[3449,]    6
[3450,]    5
[3451,]    7
[3452,]    7
[3453,]    5
[3454,]    7
[3455,]    6
[3456,]    5
[3457,]    6
[3458,]    6
[3459,]    5
[3460,]    4
[3461,]    4
[3462,]    5
[3463,]    5
[3464,]    1
[3465,]    2
[3466,]    2
[3467,]    6
[3468,]    2
[3469,]    3
[3470,]    4
[3471,]    5
[3472,]    5
[3473,]    6
[3474,]    6
[3475,]    6
[3476,]    4
[3477,]    8
[3478,]    6
[3479,]    6
[3480,]    7
[3481,]    6
[3482,]    6
[3483,]    9
[3484,]    8
[3485,]    4
[3486,]    8
[3487,]    4
[3488,]    5
[3489,]    9
[3490,]    8
[3491,]    9
[3492,]    4
[3493,]    4
[3494,]    5
[3495,]    6
[3496,]    6
[3497,]    4
[3498,]    5
[3499,]    5
[3500,]    6
[3501,]    6
[3502,]    6
[3503,]    7
[3504,]    9
[3505,]    8
[3506,]    5
[3507,]    7
[3508,]    7
[3509,]    8
[3510,]    6
[3511,]    9
[3512,]    4
[3513,]    8
[3514,]    8
[3515,]    5
[3516,]    7
[3517,]    4
[3518,]    3
[3519,]    4
[3520,]    6
[3521,]    5
[3522,]    5
[3523,]    8
[3524,]    8
[3525,]    7
[3526,]    7
[3527,]    3
[3528,]    3
[3529,]    5
[3530,]    6
[3531,]    7
[3532,]    4
[3533,]    4
[3534,]    5
[3535,]    4
[3536,]    8
[3537,]    9
[3538,]    4
[3539,]    4
[3540,]    7
[3541,]    6
[3542,]    5
[3543,]    3
[3544,]    5
[3545,]    4
[3546,]    3
[3547,]    4
[3548,]    8
[3549,]    8
[3550,]    7
[3551,]    7
[3552,]    7
[3553,]    6
[3554,]    8
[3555,]    4
[3556,]    4
[3557,]    4
[3558,]    4
[3559,]    8
[3560,]    9
[3561,]    7
[3562,]    6
[3563,]    9
[3564,]    7
[3565,]    9
[3566,]    6
[3567,]    7
[3568,]    7
[3569,]    7
[3570,]    8
[3571,]    3
[3572,]    6
[3573,]    8
[3574,]    9
[3575,]    5
[3576,]    7
[3577,]    5
[3578,]    5
[3579,]    7
[3580,]    7
[3581,]    8
[3582,]    6
[3583,]    5
[3584,]    4
[3585,]    5
[3586,]    7
[3587,]    5
[3588,]    3
[3589,]    6
[3590,]    4
[3591,]    5
[3592,]    6
[3593,]    6
[3594,]    5
[3595,]    4
[3596,]    3
[3597,]    9
[3598,]    6
[3599,]    8
[3600,]    4
[3601,]    8
[3602,]    5
[3603,]    8
[3604,]    6
[3605,]    6
[3606,]    5
[3607,]    6
[3608,]    9
[3609,]    7
[3610,]    6
[3611,]    2
[3612,]    5
[3613,]    3
[3614,]    2
[3615,]    6
[3616,]    3
[3617,]    5
[3618,]    2
[3619,]    5
[3620,]    5
[3621,]    7
[3622,]    7
[3623,]    6
[3624,]    5
[3625,]    9
[3626,]    7
[3627,]    9
[3628,]    8
[3629,]    6
[3630,]    5
[3631,]    3
[3632,]    6
[3633,]    7
[3634,]    7
[3635,]    2
[3636,]    5
[3637,]    5
[3638,]    7
[3639,]    8
[3640,]    9
[3641,]    3
[3642,]    2
[3643,]    2
[3644,]    5
[3645,]    6
[3646,]    6
[3647,]    7
[3648,]    5
[3649,]    5
[3650,]    6
[3651,]    6
[3652,]    9
[3653,]    6
[3654,]    9
[3655,]    9
[3656,]    8
[3657,]    7
[3658,]    7
[3659,]    5
[3660,]    7
[3661,]    8
[3662,]    5
[3663,]    8
[3664,]    4
[3665,]    3
[3666,]    3
[3667,]    9
[3668,]    8
[3669,]    6
[3670,]    7
[3671,]    7
[3672,]    7
[3673,]    1
[3674,]    8
[3675,]    6
[3676,]    7
[3677,]    6
[3678,]    7
[3679,]    4
[3680,]    0
[3681,]    8
[3682,]    6
[3683,]    8
[3684,]    6
[3685,]    5
[3686,]    9
[3687,]    8
[3688,]    4
[3689,]    4
[3690,]    7
[3691,]    8
[3692,]    6
[3693,]    8
[3694,]    6
[3695,]    8
[3696,]    6
[3697,]    4
[3698,]    4
[3699,]    5
[3700,]    5
[3701,]    7
[3702,]    7
[3703,]    2
[3704,]    3
[3705,]    5
[3706,]    5
[3707,]    8
[3708,]    5
[3709,]    8
[3710,]    6
[3711,]    5
[3712,]    8
[3713,]    5
[3714,]    6
[3715,]    7
[3716,]    6
[3717,]    8
[3718,]    8
[3719,]    7
[3720,]    8
[3721,]    6
[3722,]    3
[3723,]    6
[3724,]    5
[3725,]    8
[3726,]    7
[3727,]    6
[3728,]    5
[3729,]    4
[3730,]    8
[3731,]    6
[3732,]    8
[3733,]    6
[3734,]    6
[3735,]    7
[3736,]    2
[3737,]    5
[3738,]    7
[3739,]    6
[3740,]    7
[3741,]    3
[3742,]    4
[3743,]    3
[3744,]    3
[3745,]    2
[3746,]    5
[3747,]    8
[3748,]    5
[3749,]    4
[3750,]    8
[3751,]    7
[3752,]    6
[3753,]    6
[3754,]    9
[3755,]    8
[3756,]    7
[3757,]    6
[3758,]    8
[3759,]    9
[3760,]    6
[3761,]    3
[3762,]    7
[3763,]    7
[3764,]    5
[3765,]    5
[3766,]    8
[3767,]    7
[3768,]    9
[3769,]    5
[3770,]    6
[3771,]    6
[3772,]    8
[3773,]    6
[3774,]    7
[3775,]    9
[3776,]    8
[3777,]    7
[3778,]    6
[3779,]    8
[3780,]    7
[3781,]    5
[3782,]    5
[3783,]    5
[3784,]    4
[3785,]    4
[3786,]    2
[3787,]    5
[3788,]    4
[3789,]    5
[3790,]    8
[3791,]    4
[3792,]    4
[3793,]    4
[3794,]    6
[3795,]    8
[3796,]    6
[3797,]    8
[3798,]    7
[3799,]    4
[3800,]    7
[3801,]    7
[3802,]    9
[3803,]    7
[3804,]    6
[3805,]    2
[3806,]    4
[3807,]    2
[3808,]    1
[3809,]    4
[3810,]    2
[3811,]    1
[3812,]    6
[3813,]    7
[3814,]    7
[3815,]    7
[3816,]    5
[3817,]    6
[3818,]    7
[3819,]    5
[3820,]    4
[3821,]    9
[3822,]    8
[3823,]    6
[3824,]    9
[3825,]    8
[3826,]    3
[3827,]    3
[3828,]    6
[3829,]    4
[3830,]    5
[3831,]    5
[3832,]    3
[3833,]    5
[3834,]    7
[3835,]    4
[3836,]    7
[3837,]    7
[3838,]    8
[3839,]    6
[3840,]    6
[3841,]    6
[3842,]    3
[3843,]    4
[3844,]    7
[3845,]    8
[3846,]    6
[3847,]    6
[3848,]    6
[3849,]    4
[3850,]    5
[3851,]    6
[3852,]    4
[3853,]    3
[3854,]    4
[3855,]    9
[3856,]    8
[3857,]    7
[3858,]    6
[3859,]    8
[3860,]    6
[3861,]    7
[3862,]    5
[3863,]    5
[3864,]    7
[3865,]    9
[3866,]    8
[3867,]    7
[3868,]    6
[3869,]    6
[3870,]    7
[3871,]    5
[3872,]    7
[3873,]    5
[3874,]    9
[3875,]    8
[3876,]    7
[3877,]    4
[3878,]    9
[3879,]    5
[3880,]    9
[3881,]    8
[3882,]    5
[3883,]    3
[3884,]    3
[3885,]    5
[3886,]    7
[3887,]    9
[3888,]    9
[3889,]    7
[3890,]    6
[3891,]    8
[3892,]    7
[3893,]    8
[3894,]    8
[3895,]    7
[3896,]    3
[3897,]    3
[3898,]    2
[3899,]    4
[3900,]    8
[3901,]    7
[3902,]    8
[3903,]    7
[3904,]    6
[3905,]    8
[3906,]    8
[3907,]    5
[3908,]    7
[3909,]    5
[3910,]    6
[3911,]    4
[3912,]    9
[3913,]    7
[3914,]    6
[3915,]    8
[3916,]    6
[3917,]    3
[3918,]    6
[3919,]    7
[3920,]    6
[3921,]    5
[3922,]    5
[3923,]    8
[3924,]    7
[3925,]    8
[3926,]    7
[3927,]    5
[3928,]    7
[3929,]    6
[3930,]    6
[3931,]    5
[3932,]    5
[3933,]    3
[3934,]    5
[3935,]    9
[3936,]    9
[3937,]    4
[3938,]    9
[3939,]    4
[3940,]    3
[3941,]    6
[3942,]    2
[3943,]    2
[3944,]    5
[3945,]    8
[3946,]    6
[3947,]    7
[3948,]    6
[3949,]    6
[3950,]    9
[3951,]    6
[3952,]    8
[3953,]    5
[3954,]    5
[3955,]    7
[3956,]    5
[3957,]    8
[3958,]    3
[3959,]    4
[3960,]    4
[3961,]    5
[3962,]    5
[3963,]    8
[3964,]    6
[3965,]    5
[3966,]    7
[3967,]    4
[3968,]    6
[3969,]    4
[3970,]    5
[3971,]    8
[3972,]    4
[3973,]    5
[3974,]    5
[3975,]    7
[3976,]    8
[3977,]    3
[3978,]    4
[3979,]    2
[3980,]    5
[3981,]    2
[3982,]    3
[3983,]    7
[3984,]    8
[3985,]    3
[3986,]    6
[3987,]    5
[3988,]    6
[3989,]    7
[3990,]    8
[3991,]    5
[3992,]    5
[3993,]    3
[3994,]    4
[3995,]    5
[3996,]    8
[3997,]    9
[3998,]    8
[3999,]    9
[4000,]    7
data.frame(obs = ppd) %>% 
  ggplot(aes(x=obs)) +
  geom_histogram() +
  labs(x="Observed water (out of 9)")

Sampling parameters from the prior

Oops, we’ve jumped ahead of ourselves! Best practice is to simulate values from your prior first and check to see if those priors are reasonable.

m1p <-
  brm(data = list(w = 6),                            
      family = binomial(link = "identity"),          
      w | trials(9) ~ 0 + Intercept,                 
      prior(beta(1, 1), class = b, lb = 0, ub = 1),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,          
      sample_prior = "only")
samples_from_prior = as_draws_df(m1p)
samples_from_prior %>% 
  ggplot(aes(x=b_Intercept)) +
  geom_density(fill = "grey", color = "white") +
  labs(x="Proportion water", title="Prior")

Simulating observations from the prior

We may also want the PRIOR PREDICTIVE DISTRIBUTION which is the expected observiations given our prior.

prior_pd = posterior_predict(m1p)
data.frame(obs = prior_pd) %>% 
  ggplot(aes(x=obs)) +
  geom_histogram() +
  labs(x="Observed water (out of 9)")

Simulating from your priors – prior predictive simulation – is an essential part of modeling. This allows you to see what your choices imply about the data. You’ll be able to diagnose bad choices.

an aside about learning in R

At this point in the course, I’m going to start throwing a lot of code at you. Do I expect you to memorize this code? Of course not.

Do you need to understand every single thing that’s happening in the code? Nope.

But, you’ll learn a lot by taking the time to figure out what’s happening in a code chunk. Class time will frequently include exercises where I ask you to adapt code I’ve shared in the slides to a new dataset or to answer a new problem. When doing so, go back through the old code and figure out what’s going on. Run the code one line at a time. Always observe the output and take some time to look at the object that was created or modified. Here are some functions that will be extremely useful:

str() # what kind of object is this? what is its structure?
dim() # what are the dimensions (rows/columns) of this object
head() # give me the first bit of this object
str(prior_pd)
 int [1:4000, 1] 5 1 0 0 7 9 0 5 5 3 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : NULL
dim(prior_pd)
[1] 4000    1
head(prior_pd)
     [,1]
[1,]    5
[2,]    1
[3,]    0
[4,]    0
[5,]    7
[6,]    9

Continous outcomes

The globe tossing example is cute and easy to work with, but let’s move towards the kinds of variables we more frequently work with. Let’s create a model for some outcome, \(y\) that is continuous.

\[\begin{align*} y_i &\sim \text{Normal}(\mu, \sigma) \\ \mu &\sim \text{Normal}(0, 10) \\ \sigma &\sim \text{Uniform}(0, 5) \end{align*}\]

set.seed(9)
y = rnorm(n = 31, mean = 4, sd = .5)
m2 = brm(
  data = list(y=y),
  family = gaussian,
  y ~ 1,
  prior = c(prior( normal(0,10), class=Intercept),
            prior( uniform(0,5), class=sigma)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  file = here("files/models/m21.2")
)

An example: weight and height

Using the Howell data (don’t load the rethinking package because it can interfere with brms).

data("Howell1", package = "rethinking")
d <- Howell1
str(d)
'data.frame':   544 obs. of  4 variables:
 $ height: num  152 140 137 157 145 ...
 $ weight: num  47.8 36.5 31.9 53 41.3 ...
 $ age   : num  63 63 65 41 51 35 32 27 19 54 ...
 $ male  : int  1 0 0 1 0 1 0 1 0 1 ...
library(measurements)
d$height <- conv_unit(d$height, from = "cm", to = "feet")
d$weight <- conv_unit(d$weight, from = "kg", to = "lbs")
rethinking::precis(d)
             mean         sd      5.5%    94.5%      histogram
height  4.5362072  0.9055921  2.661042   5.4375      ▁▁▁▁▂▂▇▇▁
weight 78.5079631 32.4502290 20.636856 120.1583 ▁▂▃▂▂▁▁▃▅▇▇▃▂▁
age    29.3443934 20.7468882  1.000000  66.1350      ▇▅▅▃▅▂▂▁▁
male    0.4724265  0.4996986  0.000000   1.0000     ▇▁▁▁▁▁▁▁▁▇
d2 <- d[ d$age >= 18, ]

exercise

Write a mathematical model for the weights in this data set. (Don’t worry about predicting from other variables yet.)

solution

\[\begin{align*} w &\sim \text{Normal}(\mu, \sigma) \\ \mu &\sim \text{Normal}(130, 20) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]

\[\begin{align*} w &\sim \text{Normal}(\mu, \sigma) \\ \mu &\sim \text{Normal}(130, 20) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]

exercise

Simulate from your priors (parameters values and prior predictive values).

solution

Sample from your priors:

m3p = brm(
  data = d2,
  family = gaussian,
  weight ~ 1,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( uniform(0,25), class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  sample_prior = "only")

Sampling parameter estimates for the Intercept.

pairs(m3p)

This is a different (shorter) way to plot your posterior. Good things: it automatically includes the scatterplot so you can see the implications of how these parameters correlate. Bad things: not customizable and not useable when you have a lot of parameters.

Simulate values of weight.

prior_pd = posterior_predict(m3p)
dim(prior_pd)
[1] 4000  352
as.data.frame(prior_pd) %>% 
  pivot_longer(everything()) %>% 
  ggplot(aes(x=value)) +
  geom_histogram() +
  labs(x="Expected observed weights (based on prior)")

Another shorter way:

pp_check(m3p)

Fit the model

m3 = brm(
  data = d2,
  family = gaussian,
  weight ~ 1,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( uniform(0,25), class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  file = here("files/models/m21.3"))
posterior_summary(m3)
                Estimate  Est.Error         Q2.5        Q97.5
b_Intercept    99.221909 0.75720356    97.751074   100.737488
sigma          14.291987 0.55198322    13.254363    15.428183
Intercept      99.221909 0.75720356    97.751074   100.737488
lprior         -8.318377 0.05827127    -8.433538    -8.203915
lp__        -1441.294276 1.02695551 -1444.235605 -1440.292326
pairs(m3)
Code
posterior_predict(m3) %>% 
  as.data.frame() %>% 
  pivot_longer(everything()) %>% 
  ggplot(aes(x=value)) +
  geom_density(fill = "grey", color = "white") +
  geom_density( aes(x = weight), data=d2, inherit.aes = F) 
pp_check(m3)

Adding in a linear component

We might assume that height and weight are associated with each other. Indeed, within our sample:

plot(d2$weight ~ d2$height)

exercise

Update your mathematical model to incorporate height.

\[\begin{align*} w_i &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha + \beta h_i \\ \alpha &\sim \text{Normal}(??, ??) \\ \beta &\sim \text{Normal}(0, 25) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]

exercise

Update your mathematical model to incorporate height.

\[\begin{align*} w_i &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha + \beta (h_i - \bar{h}) \\ \alpha &\sim \text{Normal}(130, 20) \\ \beta &\sim \text{Normal}(0, 25) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]

\[\begin{align*} w_i &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha + \beta (h_i - \bar{h}) \\ \alpha &\sim \text{Normal}(130, 20) \\ \beta &\sim \text{Normal}(0, 25) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]

To update our brms code:

d2$height_c = d2$height - mean(d2$height)
  
m4p = brm(
  data = d2,
  family = gaussian,
  weight ~ 1 + height_c,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( normal(0,25),   class=b),
            prior( uniform(0,25),  class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  sample_prior = "only")
set.seed(9)
samples_from_prior = as_draws_df(m4p)
str(samples_from_prior)
draws_df [4,000 × 9] (S3: draws_df/draws/tbl_df/tbl/data.frame)
 $ b_Intercept: num [1:4000] 167 155.6 149.6 88.7 176.4 ...
 $ b_height_c : num [1:4000] -57.9 -57.2 -18.5 52 -47.3 ...
 $ sigma      : num [1:4000] 14.1 11.8 13.3 2.1 21.1 ...
 $ Intercept  : num [1:4000] 167 155.6 149.6 88.7 176.4 ...
 $ lprior     : num [1:4000] -15.7 -14.7 -12 -15.6 -15.8 ...
 $ lp__       : num [1:4000] -13.8 -12.9 -10.2 -14.9 -14.6 ...
 $ .chain     : int [1:4000] 1 1 1 1 1 1 1 1 1 1 ...
 $ .iteration : int [1:4000] 1 2 3 4 5 6 7 8 9 10 ...
 $ .draw      : int [1:4000] 1 2 3 4 5 6 7 8 9 10 ...
Code
d2 %>% 
  ggplot(aes(x=height_c, y=weight)) +
  geom_blank() +
  geom_abline( aes(intercept=b_Intercept, slope=b_height_c), 
               data=samples_from_prior[1:50, ],
               alpha=.3) +
  scale_x_continuous(name = "height(feet)", 
                     breaks=seq(4,6,by=.5)-mean(d2$height),
                     labels=seq(4,6,by=.5))

Describe in words what’s wrong with our priors.

Slope should not be negative. How can we fix this?

Could use a uniform distribution bounded by 0.

m4p = brm(
  data = d2,
  family = gaussian,
  weight ~ 1 + height_c,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( uniform(0,25),   class=b),
            prior( uniform(0,25),  class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  sample_prior = "only")

exercise

Fit the new weight model to the data.

solution

m4 = brm(
  data = d2,
  family = gaussian,
  weight ~ 1 + height_c,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( uniform(0,25),   class=b),
            prior( uniform(0,25),  class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  file=here("files/models/m21.4"))
posterior_summary(m4)
               Estimate Est.Error        Q2.5       Q97.5
b_Intercept    99.21508 0.5154139    98.17317   100.20955
b_height_c     24.74661 0.2609169    24.05150    24.99355
sigma          10.36090 0.3898531     9.65234    11.12228
Intercept      99.21508 0.5154139    98.17317   100.20955
lprior        -11.53739 0.0396881   -11.61861   -11.46176
lp__        -1332.15882 1.3936194 -1335.90639 -1330.52002

exercise

Draw lines from the posterior distribution and plot with the data.

solution

Code
set.seed(9)
samples_from_posterior = as_draws_df(m4)
d2%>% 
  ggplot(aes(x=height_c, y=weight)) +
  geom_point(size=.5) +
  geom_abline( aes(intercept=b_Intercept, slope=b_height_c), 
               data=samples_from_posterior[1:50, ],
               alpha=.3,
               color="#1c5253") +
  scale_x_continuous(name = "height(feet)", 
                     breaks=seq(4,6,by=.5)-mean(d2$height),
                     labels=seq(4,6,by=.5))

A side note: a major concern or critique of Bayesian analysis is that the subjectivity of the priors allow for nefarious behavior. “Putting our thumbs on the scale,” so to speak. But priors are quickly overwhelmed by data. Case in point:

m4e = brm(
  data = d2,
  family = gaussian,
  weight ~ 1 + height_c,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( normal(-5,5),   class=b),
            prior( uniform(0,25),  class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  file=here("files/models/m21.4e"))
posterior_summary(m4e)
                Estimate Est.Error         Q2.5       Q97.5
b_Intercept    99.197733 0.5035653    98.214092   100.15932
b_height_c     35.775818 1.8832462    32.177572    39.50070
sigma           9.529429 0.3687178     8.839379    10.27233
Intercept      99.197733 0.5035653    98.214092   100.15932
lprior        -44.172476 3.0738968   -50.456409   -38.49994
lp__        -1334.629214 1.1849503 -1337.645046 -1333.23509

You’ll only really get into trouble with uniform priors that have a boundary, if true population parameter is outside your boundary. A good rule of thumb is to avoid the uniform distribution. We’ll cover other options for priors for \(\sigma\) in future lectures, but as a preview, the exponential distribution works very well for this!